SETUP

Load packages and paths

library(grid)
library(ggplot2)
library(plyr)
library(RColorBrewer)

Color Scheme

CMI Recommendations

cmi_main_blue = "#0071b2"
cmi_grey = "#929d9e"
cmi_light_blue = "#00c4d9"
cmi_pea_green = "#b5bf00"

cmi_rich_green = "#73933d"
cmi_rich_purple = "#8e7fac"
cmi_rich_red = "#d75920"
cmi_rich_blue = "#4c87a1"
cmi_rich_aqua = "#66c7c3"
cmi_rich_orange = "#eebf42"

cmi_vibrant_yellow = "#ffd457"
cmi_vibrant_orange = "#f58025"
cmi_vibrant_green = "#78a22f"
cmi_vibrant_garnet = "#e6006f"
cmi_vibrant_purple = "#9A4d9e"
cmi_vibrant_blue = "#19398a"

cmi_site_colors = c(cmi_vibrant_blue, cmi_rich_blue, cmi_vibrant_purple, cmi_vibrant_garnet, 
    cmi_rich_red, cmi_vibrant_orange, cmi_vibrant_yellow, cmi_vibrant_green)
cmi_site_colors_ramp = colorRampPalette(cmi_site_colors)

Load data

Read in the data and then some

# setwd('/home2/zarrar/projects/qc') setwd('~/zarrar/qc')
# setwd('~/Dropbox/Research/cmi/qc')
df <- read.csv("../corr.qc/qc_filt_epi_derivatives.csv")[, -1]
nsites <- length(unique(df$site))

Percentiles

In our plots, we want to have percentile lines for each QC measure to indicate the distribution of each site relative to the whole sample

qc.measures <- colnames(df)[!(colnames(df) %in% c("uniqueid", "subid", "site", 
    "site.name", "session", "scan", "global"))]
qvals <- c(0.01, 0.05, 0.25, 0.5, 0.75, 0.95, 0.99)
qcat <- c(1, 5, 25, 50, 25, 5, 1)
qline <- c(3, 2, 5, 1, 5, 2, 3)
qsize <- c(0.4, 0.25, 0.3, 0.25, 0.3, 0.25, 0.4)
qcols <- c("grey10", "grey10", "grey10", "grey50", "grey10", "grey10", "grey10")
# qcols <- brewer.pal(8, 'Dark2')[c(1,2,3,4,3,2,1)]

Now let's get the percentiles

percentiles <- apply(subset(df, select = qc.measures), 2, quantile, qvals, na.rm = TRUE)
percentiles <- as.data.frame(cbind(percentiles, qcat, qline, qsize))
percentiles$qline <- as.factor(qline)
percentiles$qcat <- as.factor(qcat)
print(percentiles)

##     falff_50 reho_50 vmhc_50 falff_75 reho_75 vmhc_75 falff_90 reho_90
## 1%    0.3900 0.05873  0.1776   0.4347 0.07763  0.3422   0.4817  0.1051
## 5%    0.5027 0.06914  0.2355   0.5555 0.09768  0.4187   0.6112  0.1349
## 25%   0.6148 0.09503  0.3212   0.6559 0.14041  0.5167   0.7014  0.1928
## 50%   0.6422 0.11562  0.3913   0.6928 0.17240  0.5821   0.7473  0.2334
## 75%   0.6849 0.15442  0.4562   0.7343 0.22562  0.6423   0.7836  0.2975
## 95%   0.7497 0.21802  0.5579   0.7844 0.31196  0.7275   0.8191  0.4010
## 99%   0.7619 0.25384  0.6079   0.8014 0.35066  0.7767   0.8364  0.4430
##     vmhc_90 falff_fwhm reho_fwhm vmhc_fwhm falff_mean reho_mean vmhc_mean
## 1%   0.5338      2.423     5.817     5.393     0.3971   0.06685    0.2136
## 5%   0.6078      2.741     8.109     5.803     0.5098   0.08103    0.2625
## 25%  0.6859      3.379    10.299     6.256     0.6203   0.11061    0.3291
## 50%  0.7370      4.157    11.582     6.465     0.6508   0.13198    0.3796
## 75%  0.7836      4.822    12.685     6.649     0.6895   0.17173    0.4256
## 95%  0.8453      5.689    13.834     6.900     0.7492   0.23689    0.4957
## 99%  0.8755      6.137    15.088     7.081     0.7614   0.26737    0.5303
##     qcat qline qsize
## 1%     1     3  0.40
## 5%     5     2  0.25
## 25%   25     5  0.30
## 50%   50     1  0.25
## 75%   25     5  0.30
## 95%    5     2  0.25
## 99%    1     3  0.40

Measure Descriptions

Associate a detailed description with each measure

dnames <- c("fALFF", "REHO", "VMHC")
mnames <- c("Median", "75th Percentile", "90th Percentile", "FWHM (mm)", "Mean")
descs <- expand.grid(m = mnames, d = dnames)
descs <- paste(descs[, 2], descs[, 1], sep = " - ")
mdf <- data.frame(measure = qc.measures, description = descs)
print(mdf)

##       measure             description
## 1    falff_50          fALFF - Median
## 2     reho_50 fALFF - 75th Percentile
## 3     vmhc_50 fALFF - 90th Percentile
## 4    falff_75       fALFF - FWHM (mm)
## 5     reho_75            fALFF - Mean
## 6     vmhc_75           REHO - Median
## 7    falff_90  REHO - 75th Percentile
## 8     reho_90  REHO - 90th Percentile
## 9     vmhc_90        REHO - FWHM (mm)
## 10 falff_fwhm             REHO - Mean
## 11  reho_fwhm           VMHC - Median
## 12  vmhc_fwhm  VMHC - 75th Percentile
## 13 falff_mean  VMHC - 90th Percentile
## 14  reho_mean        VMHC - FWHM (mm)
## 15  vmhc_mean             VMHC - Mean

cat("CHECK ABOVE...do the columns for each row match\n")

## CHECK ABOVE...do the columns for each row match

PLOTS

More Setup

This function will add percentile lines in the background plot: ggplot object pdf: percentile data frame

compile_percentiles <- function(pdf, measure, cols = NULL) {
    ret <- lapply(1:nrow(pdf), function(i) {
        p <- pdf[i, ]
        if (!is.null(cols)) {
            plot <- geom_hline(aes_string(yintercept = measure), data = p, size = as.numeric(p$qsize), 
                linetype = as.numeric(p$qline), color = cols[i])
            # as.character(p$qcolor[1])
        } else {
            plot <- geom_hline(aes_string(yintercept = measure), data = p, size = as.numeric(p$qsize[1]), 
                linetype = as.numeric(p$qline[1]), color = "grey50")
        }
        return(plot)
    })
    return(ret)
}

Outliers

Sometimes extreme data-points can skew the plot and make it difficult to see the spread of the data.

I will avoid plotting those outlier points by setting the axis to only the points that I want. ggplot will complain but let her (poor baby).

# functions
range.outlier.iqr <- function(x, times = 3) {
    upper.limit <- quantile(x, 0.75) + times * IQR(x)
    lower.limit <- quantile(x, 0.25) - times * IQR(x)
    return(c(lower.limit, upper.limit))
}
outlier.iqr <- function(x, times = 3) {
    tmp <- range.outlier.iqr(x, times)
    lower.limit <- tmp[1]
    upper.limit <- tmp[2]
    return((x > upper.limit) | (x < lower.limit))
}
# outlier values (if any)
lst.outlier.iqr <- llply(qc.measures, function(measure) {
    ret <- subset(df, select = c("uniqueid", "subid", "site", "site.name", "session", 
        "scan", measure))
    inds <- outlier.iqr(df[[measure]])
    return(ret[inds, ])
})
names(lst.outlier.iqr) <- qc.measures
# new ranges of our plots (sans outliers)
df.range.iqr <- as.data.frame(sapply(qc.measures, function(m) {
    inds <- !outlier.iqr(df[[m]])
    range(df[[m]][inds]) * c(0.99, 1.01)
}))

Visualization of Text

A function with all the theme jazz

set_themes <- function(family = "Times", text.size.x = 14, text.size.y = 16, 
    title.size = 18) {
    family <- "sans"
    pg <- list(theme_bw(), theme(axis.title.x = element_text(family = family, 
        face = "plain", size = title.size)), theme(axis.title.y = element_text(family = family, 
        face = "plain", size = title.size, angle = 90, vjust = 0.25)), theme(axis.text.x = element_text(family = family, 
        face = "plain", size = text.size.x, vjust = 0.5, angle = 45)), theme(axis.text.y = element_text(family = family, 
        face = "plain", size = text.size.y, angle = 90)), theme(axis.ticks.length = unit(0.15, 
        "lines")), theme(axis.ticks.margin = unit(0.15, "lines")), theme(plot.margin = unit(c(0.25, 
        1, 0.25, 1), "lines")), theme(legend.position = "none"))
    return(pg)
}

QC Derivative Measures

I will be plotting a bunch of different data-sets here. First I will have all the data, then I will have it remote 3 x the IQR. So two sets of the same plots.